#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
This module provides the Mammoth and MammothControls classes in order to
master the mammoth !
The MammothControls class store all the mammoth parameters and do the needed
checks.
The Mammoth class provides methods in order to setup the calculation or the
wanted information. Each mammoth features is associated to a mammoth method.
The Mammoth class use the MammothControls class in order to do what you want.
"""
from pymatgen.io.gaussian import GaussianOutput
from orca import OrcaHessian, OrcaOutfile, OrcaEnGradfile
from .core import Molecule
from .qm2ff import QM2FF
__all__ = ["Mammoth", "MammothControls"]
# Ha.bohr^-2 into kcal.mol^-1 A^-2
HA_BOHR_m2_TO_KCAL_MOL_m1_A_m2 = 2240.874995
[docs]class Mammoth:
""" Mammoth class to master the Mammoth """
# TODO: at First there is mainly a main methods
# Next, one can imagine a specific method for the various purpose or
# type of Mammoth calculations
def __init__(self, controls):
"""
Set up the mammoth to run the desired calculations
Args:
controls (MammothControls): The MammothConrols object that contains
all the control parameters
"""
self.controls = controls
self.molecule = None
[docs] def setup_molecule(self):
""" Setup a the molecule object from a QM results"""
if self.controls.software_qm == "orca":
# read in ORCA out files
ohess = OrcaHessian(self.controls.filename + ".hess")
oout = OrcaOutfile(self.controls.filename + ".out")
ograd = OrcaEnGradfile(self.controls.filename+ ".engrad")
species = ohess.molecule.species
coords = ohess.get_mol_ang().cart_coords
gradient = ograd.gradient
# select atomic charges
if self.controls.charges_partition == "mulliken":
charges = oout.mul_charges
elif self.controls.charges_partition == "loewdin":
charges = oout.loe_charges
# symmetrize Hessian
if self.controls.symmetrize_hessian:
# TODO: Comment about the -1 of the hessian matrix ?
hessian = -1 * ohess.symmetrize_hessian()
else:
hessian = ohess.hessian
# convert to kcal.mol-1.A-2
hessian *= HA_BOHR_m2_TO_KCAL_MOL_m1_A_m2
bond_orders = oout.bond_orders
elif self.controls.software_qm == "gaussian":
gout = GaussianOutput(self.controls.filename + ".out")
species = gout.final_structure.species
coords = gout.final_structure.cart_coords
hessian = -1 * gout.hessian
# convert to kcal.mol-1.A-2
hessian *= HA_BOHR_m2_TO_KCAL_MOL_m1_A_m2
charges = [data[1] for data in gout.Mulliken_charges.values()]
# if bond orders exist, consider only those larger than .1 (as done
# in orca)
if gout.bond_orders:
bond_orders = {k: v for k, v in gout.bond_orders.items() if v > .1}
else:
print("WARNING: No bond orders in gaussian output file %s" %
self.controls.filename)
bond_orders = {}
# set up the molecule object
self.molecule = Molecule(species=species, coords=coords, hessian=hessian,
bond_orders=bond_orders, at_charges=charges, grad=gradient)
[docs] def main(self):
""" Run the FF calculation """
print(self.controls)
# set up the molecule
self.setup_molecule()
run = QM2FF(self.molecule, **self.controls.get_qm2ff_params())
forcefield = run.get_forcefield()
print(forcefield)
if self.controls.software_md == "lammps":
forcefield.export_lammps_data(self.controls.filename + ".data")
[docs]class MammothControls(dict):
""" A convenient container of all Mammoth control paramters """
# list of authorized parameters with default values
# the rule is that all keys and values are lowercase
parameters = {
"filename": "",
# QM parameters
"software_md": "lammps",
"software_qm": "orca",
"charges_partition": "loewdin",
# FF parameters
"bond_search": "hessian",
"cutoffb": 2.5,
"symmetrize_hessian": "true",
"dihedral_type": "harmonic",
"improper_type": "all",
"cutoff_sp2": 0.15,
# MD parameters
"temperature": 100.,
"timestep": 0.25,
"total_time": 1
}
# parameters type: use to convert parameter values in the given type
parameters_type = {
"filename": str,
# QM parameters
"software_md": str,
"software_qm": str,
"charges_partition": str,
# FF parameters
"bond_search": str,
"cutoffb": float,
"symmetrize_hessian": bool,
"dihedral_type": str,
"improper_type": str,
"cutoff_sp2": float,
# MD parameters
"temperature": float,
"timestep": float,
"total_time": float
}
# parameters authorized values if relevant
parameters_values = {
"software_md": ("lammps", "gromacs"),
"software_qm": ("orca", "gaussian"),
"charges_partition": ("mulliken", "loewdin"), # "mk"),
"bond_search": ("hessian", "bond_order", "hybrid", "connectivity"),
"symmetrize_hessian": ("true", "false"),
"dihedral_type": ("harmonic", "parabolic"),
"improper_type": ("all", "sp2")
}
def __init__(self, *args, **kwargs):
# first, set up default values
self.update(self.parameters)
# now update with given Args
self.update(*args, **kwargs)
# check if filename is empty
if self["filename"] == "":
raise ValueError("You must specify a filename.")
[docs] def update(self, *args, **kwargs):
# update parameters all checks are done
# following the __setitem__ method
for k, v in dict(*args, **kwargs).items():
self[k] = v
[docs] @classmethod
def from_file(cls, inputfile="input.mmt", encoding="utf8"):
""" read in inputfile and load all parameters in a new instance """
params = {}
with open(inputfile, "r", encoding=encoding) as f:
for line in f:
if line[0] in {"#", "!", "&"} or line.strip() == "":
# skip blanc line or line starting by a special character
continue
else:
# read key = value line
k, v = line.split()[0:3:2]
# all parameters are stored in lowercase
params[k.lower()] = v.lower()
return cls(**params)
[docs] @classmethod
def default_parameters(cls):
"""
All default parameters as a dict that can be used to set up a
MammothConrols object. The filename is empty.
"""
return cls.parameters.copy()
[docs] def get_qm2ff_params(self):
""" return a dict with the needed parameters for the QM2FF class """
params = ["cutoffb", "bond_search", "dihedral_type", "improper_type",
"cutoff_sp2"]
return {k: self[k] for k in params}
[docs] def get_md_params(self):
""" return a dict with the needed parameters for the MD run """
params = ["temperature", "timestep", "total_time"]
return {k: self[k] for k in params}
def __setitem__(self, name, value):
""" set a parameters following dict syntax """
add = False
name = name.lower()
if isinstance(value, str):
value = value.lower()
# check if name is an available parameters
if name in self.parameters:
# check if an authorized values list is provided
if name in self.parameters_values:
if value in self.parameters_values[name]:
add = True
else:
m = "Value '%s' is not an authorized value" % value
m += " for parameter '%s'\n" % name
m += "authorized values are:\n"
m += "\n".join(str(val) for val in self.parameters_values[name])
raise ValueError(m)
else:
add = True
else:
raise KeyError("Parameters '%s' does not exist." % name)
if add:
instance = self.parameters_type[name]
try:
if instance == bool:
bool_value = True if value.lower() == "true" else False
super().__setitem__(name, bool_value)
else:
super().__setitem__(name, instance(value))
except ValueError:
raise ValueError("Could not convert parameter '%s' to %s" % (name, instance))
def __getattr__(self, name):
""" get paramters as an attribute of the class """
return self.get(name)
def __str__(self):
""" Return a string with a table of all paramters """
line = 46 * "-" + "\n"
line += "Mammoth parameters".center(46) + "\n"
line += 46 * "-" + "\n"
line += "Parameters".center(20) + " " + "values".center(20) + "\n"
line += 46 * "-" + "\n"
for k, v in self.items():
line += k.rjust(20) + " = " + str(v).ljust(20) + "\n"
line += 46 * "-" + "\n"
return line
def __repr__(self):
""" Return a string representation of the instance """
line = self.__class__.__name__
for i, (k, v) in enumerate(self.items()):
if i == 0:
if isinstance(v, str):
line += "(%s='%s'" % (k, str(v))
else:
line += "(%s=%s" % (k, str(v))
else:
if isinstance(v, str):
line += ", %s='%s'" % (k, str(v))
else:
line += ", %s=%s" % (k, str(v))
line += ")"
return line